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INTRODUCTION 


Flow past a body is, in general, specified by a variety of parameters such as 
thickness, angle of attack, camber, Mach number etc. A particular flow is, therefore, 

characterized by a single point in the corresponding parameter space. Conversely, the 

numerical calculation of a particular flow field yields information at just one point of 
the parameter space. However, the nature of a continuous range of nearby flow fields 
is of fundamental significance in the design and performance of aircraft. To treat this 
generally, one can consider the variational equations (which are linear) obtained by 
differentiating the exact equations with respect to each of the relevant parameters. The 
resulting matrix of derivatives of flow quantities is referred to as the Jacobi matrix. 

The subsequent procedure is, in principle, straightforward. One integrates the 
nonlinear governing equations - which results in the determination of just one point in 
parameter space - and simultaneously the variational equations governing the Jacobi 
matrix. The last is then used to describe the neighborhood of the already determined 
point of the parameter space. A method is presented herein which allows efficient 

generation of solutions in the neighborhood of a base solution. Since the variational 
equations are linear, the additional computational time required for their integration is 

modest. 

We have applied the Jacobi matrix technique to the direct calculation of 
inviscid supersonic flow about 

o two dimensional airfoils of varying thickness, angle of attack and camber 
o axisymmetric bodies of varying thickness and taper 
and the design (inverse) calculation of inviscid supersonic flow past 

o airfoils described by a given family of pressure distributions 
o axisymmetric bodies described by a given family of pressure distributions. 
Also we applied the method to subsonic potential flow about two dimensional 
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airfoils by modifying Jameson’s FL036. 

Results of our calculations show that the Jacobi method allows for the 
efficient and accurate generation of parametric solutions in the neighborhood of a 
known solution. In general terms, we consider a system of nonlinear partial differential 
equations 

F (u(x;e),x;€) = 0 (LI) 

in the flow variables u, dependent variables x, and parameters e_. For purposes of 
exposition we regard u(x;ej as known and seek the solution at a neighboring point in 
parameter space. The parametrically differentiated dependent variables are governed by 
the equations obtained by differentiating (1), viz. 


f Fj(u *;€> fp fij + |Ei = 0 


d€ 


9Uj 0€ k 3e k 


( 1 . 2 ) 


The, in general, non-square matrix f^-i is known as the Jacobi matrix and the above 

9£ k 

procedure provides a linear system of equations governing the Jacobi matrix. The term 
9Fj/0Uj actually represents an operator, the details of which are best left to the 
individual cases. If u°= ufeig) represents a known solution of the flow then any 
neighboring flow at some fixed point x is determined by 

0U° 


u(x;e)*u°+ — (i^ i-o) 


(1.3) 


In what follows we will be somewhat loose in not distinguishing between the two sides 
of (1.3). A basic difficulty with what has been just said, in particular to the use of 
(1.3), is the fact that the conditions on the problem occur at locations which vary with 
€_. Specifically, both the boundary locations (and shock locations) may vary with changes 
in the parameters e.. We first present a method that avoids the difficulties implicit in 
such spatial variations with £, and later treat directly the formulation implicit in 
(1.1-3). 


- 2 - 


Chapter I 

The Jacobi Matrix Technique and 
It’s Application to Two-Dimensional 
Supersonic Flow 


The White Rabbit put on his spectacles. "Where shall I begin, please your 
Majesty?" he asked. 

"Begin at the beginning," the King said, very gravely, "and go on till you come to 
the end: then stop." 

- Alice’s Adventures in Wonderland 
Lewis Carroll 
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APPLICATION TO 2D SUPERSONIC FLOW 


To illustrate this method we consider steady, inviscid, supersonic flows past two 
dimensional airfoils. For this purpose and in order to be specific, consider a family 
of profiles depending on three parameters (thickness, camber, and angle of attack). For 
completeness, we summarize the methods used in solving such flows [1], [2]. The 
equations are written in characteristic form as follows: 


S 3=° 


(i.i) 

(9 + P(M) ) a = 

sin2ft s - ( a) 
27 v 

(1.2) 

(0 - P(g) ) 3 = 

x 3 

(1 - tan0 tan jz) ~ 
A a 

(1.3) 


Here the coordinates (a,0) correspond to the streamlines, a = constant, and the C + 
characteristics, 0 = constant (Figure 1). 0 is the flow deflection angle, g is the Mach 
angle and s is the entropy. P(g) is the Prandtl function given by 

P(g) = X^ tan' 1 (X^tang) - g , X = (7+1)/ (7-1) . (1.4) 

An advantage to solving the above characteristic form of the equations is that it 
generates a body fit, shock fit coordinate system. We mention in passing that since the 
equations are exact, they are valid in the hypersonic flow regime so long as such real 
gas effects as disassociation and ionization can be ignored. 

The physical coordinates x,y satisfy the relations [1], [2] 

y a = x « ^ tan ( e+ *0. y o = x e tan 0 (1.5) 

The transformation to (a,0) coordinates leaves open two arbitrary functions and these 
are fixed so that the shock is along a = 0 and the airfoil is positioned along the line 
a = 0 (Figure 2). Appropriate boundary conditions at the body are 
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( 1 . 6 ) 


x((0,8) =8, y(0,8) = f(B,e), 6(0,8) = tan^B, € )) . 

The Ranldne-Hugoniot conditions govern the jumps in 6, f l and s at the shock. 
Written in terms of the shock angle n, they are given by 

1 (M 2 - l)t an 2 n - 1 

tanG = — - — - 

tanI? (1+ 2±1 m 2 ) + (1+^1 M 2 ) tan 2 h (1.7) 

2 2 


sin 2 /r = 


0.2(l+-w )(1 + i w) 
6 6 


(1+w) (1+0.2M ) - (l+-w) (1+ iw) 
6 6 


( 1 . 8 ) 


where 


s = 2.5 In (1 + -w) + 3.5 In (1+lw) - 3.5 In (1+w) 
6 6 


w = M 2 sin 2 n - l 


(1.9) 


( 1 . 10 ) 


The shock angle is related to the coordinates as follows 

y« + ye 


tanb = dy / dx I 


shock x a + x 8 ^ shock 


( 1 . 11 ) 


In the above we have assumed a perfect gas with constant specific heats and 
hence that 

p =p? exp [( 7 - l)s] (1.12) 

It should be noted that this formulation eliminates the difficulty mentioned in the 
Introduction. Namely, by using the (a,8) - coordinate system, a quantity such as 

3p 

gM«,8 ; D 

signifies the variation with e. at fixed a and 8. In particular it gives the variation of 
pressure say fixed at the body, a = 0, or at the shock, a = 8. This makes the 
integration of the differential equations significantly simpler. 
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VARIATIONAL EQUATIONS 


We are interested in solutions to these equations at points near a known solution. 
To pursue this we differentiate all of the above equations with respect to a typical 
parameter of interest. In keeping with the remarks at the close of the previous 

section, we emphasize that differentiation is with respect to ± with a and 6 held 
fixed. 


The mechanics of the differentiation are straightforward but tedious. We represent 
by capitalized variables the differentiated variables; 


0 


30 3s 

» S = — 

de de 


X 


3x 3y 3 /j l 

> Y = __ > Y a — * 

de de de 


(1.13) 


when (1.4), (1.5), (1.6), and (1.8) are parametrically differentiated we obtain, 

s 0 = 0 


(1.14) 


©a + PVW S0cC ° S2 M 1 ▼+P (t 00* a - Si n2M S a = 0 

7 2 


(1.15) 


e 3 + 


x 

- — sec 2 0 tan 11 


0 - (1 - tan 0 tan/t) X B e _ 


a 


Xo0 t X J3 

PfjXtf) % + -2-2- sec 2 Mtan6]'f+ (1- tan0 tanfi) — - (Xg X a ) (1.16) 

x a x oc x a 


Y a = X a tan( 6+ ,*) + X a sec 2 (0 + jt)(e + V) 


(1.17) 


Yg = Xg tan 0 + Xg sec 2 0 (1-18) 

It should be noted that we have dropped the specification that e be a vector. This 

has been done for ease of exposition. This can be done without loss of generality. 
Variation with respect to each parameter can be treated separately, since only first 
order variations are being considered. 


- 6 “ 




At the shock the parametrically differentiated equations are 


dr? cos 


de ( x a +x e)' 


— 2 [ Y a( x a +x 0) + Y g( x a +x - X X B ] 


(1.19) 


17.5 3.5 3.5 

6 +7w 6 +w 1+w 


dw 

de 


( 1 . 20 ) 


0. 2A [1* 2 w) - 02 (i.2»)(i. 1w)(02M 2 - 1. 1.) 


Y = 


3 18 J dw 


A 2 sin2/x 


— ( 1 . 21 ) 

de 


dw d n -.7 1 

where ^7 = M 2 sin 2n » A = l+0.2M 2 (l+w) - (1+^w) (1+ ±w) 


dr? - f 1 

= COS^Ot- - y 

de L sin 7 ? 


d 0 dr? 
de 


(M 2 -l)tan 2 r? -1 


(1+ 2±i M 2 )+ (1+ M ^) tan 2 r? 


[(1+^±1m 2 ) + (1+^ M 2 )tan 2 r?] 
2 2 


1 [(M 2 - 1) sec r?(l+ 2+1 m 2 ) 


( 1 . 22 ) 


+ ( 1 + M 2 )tan 2 r? ] - [ (M 2 -l)tan 2 r?- lj (1+ 221 M 2 )sec 2 r?J | 


In the actual integration (1.21-1.25) are applied at the shock 
a = 3 

At the body a = 0 the appropriate equations are 

0/3 0f o(B, e ) , 

X = 0, Y = f ,0= — sr cos 2 0 

0 e 0 e 


(1.23) 


(1.24) 


In writing (1.24) we revert to the general case in which many parameters are 
being considered. At this point we can simultaneously numerically integrate the 
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non-linear system and the variational equations. The calculation of the base flow is 
second order accurate [2]. The calculation of the new flow is first order in space, 
second order in the parameters of interest. The calculation of the two flows is 
interleaved in that after the flow along 0 = constant is computed by the base code, 
the parametric code then calculates the exact derivatives in order to obtain the 
variational flow. 

RESULTS 

As we have already mentioned the method applies generally to many independently 
varying parameters. As a typical use of the variational quantities we use Taylor’s 
theorem to consider the change in pressure, 

Pnew * Phase + (dp /9e i0 ) ( e p € j0 ) (1-25) 


where € ; represents the various parameters with the differential coefficients calculated 
holding a and 0 fixed and the zero subscript denotes a reference or base calculation. 

For example, if just thickness is considered and denoted by say e, then at the 
body 

p * p + &- ] ( e ' e n ) = p + rr 

new base *-06 -'a=O,0 0 base 9e '*x,y=f 


rd P 

r~ 

CO 

V- 

*-0y 

^ x, y= f ^ 


0 > ( 1 - 26 ) 

The second form exhibits the result obtained if variation in the physical plane is 
considered. 

In the numerical calculations that are discussed, we have taken for f a family of 
shapes given by 

y = 2ex(l-x) - x tan A + lOxe (x- 1) (x - i ) c (1.27) 


Thus, e is the thickness ratio based on chord, A the mean chord angle of attack, and 
c a scaling factor for the camber (shape) function. Figure 3 shows the effect of 
changing just the thickness (A,c = 0). Here we have plotted the pressure distribution 
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on the upper surface and the airfoil which, for aesthetic reasons, has the lower surface 
plotted as a reflection of the upper surface. Note that the method gives good 
agreement with the exact solution even when the new thickness is fifty percent greater 
than the base thickness. 

More generally we consider variations in all three parameters. Thus, the pressure 
relation at the body is 


p ** p + 

new base 





(1.28) 


Figures 4 through 6 show the effect of changing various combinations of thickness, 
angle of attack, and camber. Here we see that, although the airfoil configurations are 
markedly different, there is very good agreement between the parametrically generated 
pressure distribution and the exact pressure distribution for the new airfoil. 


INVERSE CASE 

The method which has been presented also works as weli on the inverse or design 
problem where the pressure on the body is known, but the shape of the body is to be 
determined. Using the Bernoulli equation and the perfect gas law one may show [1] 


“ ■ sin ' 1 (pr 


7 - 1 

exp (-y (s + lnyp)) 


7-1 2 7-1 

1 + — M - exp( — (s +1 nyp)) 


f] 


(1.29) 


This when differentiated, yields 


Y - 


where 


(7-D 2 


O + ^M 2 


2:1 2 


2 M )exp(w) s = ( 7- l) z ( 1 + 2 M ) exp(w) 


^ [l+2^M 2 - exp(w) ] 2 sin 2 g [l+Z^M - exp(w)] 2 sin 2 g 


d(lnp) 

de 

(1.30) 


7-1 

w = — y (s+ln p) 
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Obviously at the airfoil we can no longer use (1.6) since we are hoping to 
determine the shape of the airfoil. Instead we must use (1.5) and (1.29). Therefore, 
the parametrically differentiated equations (1.24) must be replaced by (1.17), (1.18) and 
(1.30). The integration may now proceed as in the direct case [2]. The results of the 
variational calculations are presented in Figures 7 through 9. Notice that even for a 
20% change in the logarithm of the pressure (corresponding to a 50% increase in 
thickness), the difference between the exact airfoil shape and the computed shape is 
less than 4%. 
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FIGURE CAPTIONS 


Fig. 1. 

Fig. 2. 
Fig. 3. 

Fig. 4. 

Fig. 5. 

Fig. 6. 

Fig. 7. 
Fig. 8. 
Fig. 9. 


Body, C+ characteristics, streamlines and C characteristics 
(dashed) in physical (x,y) plane, from [1]. 

Body, C+ characteristics and streamlines in (a,/3) plane. 

Pressure distribution on 10% and 15% thick airfoils at M = 2 and 

10% and 15% airfoils. 

Pressure distribution on 10% thick airfoil, uncambered at 0 degree 

angle of attack and cambered (c=0.1) at 5 degrees angle of attack 

at M = 2 along with respective bodies. 

Pressure distribution on 10% thick airfoil, uncambered at 0 degree 

angle of attack and canbered (c=0.2) at 5 degrees angle of attack 

at M = 2 along with respective bodies. 

Pressure distribution on 10% thick airfoil, uncambered at 0 degree 

angle of attack and cambered (c=0.2) at 10 degrees angle of 
attack at M = 2 along with respective bodies. 

Inverse case: pressure distribution on 10% and 12% airfoils, M = 2, 

along with generated bodies. Dashed airfoil is computed shape. 

Inverse case: pressure distribution on 10% and 15% airfoils, M = 2, 

along with generated bodies. Dashed airfoil is computed shape. 

Inverse case: pressure distribution on 10% and 12% airfoils, M = 4, 

along with generated bodies. Dashed airfoil is computed shape. 
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PHYSICAL PLANE 



Figure 1 


COMPUTATIONAL PLANE 
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Chapter II 

The Application of the Jacobi Matrix Technique 
to Axisymmetric Supersonic Flow 


"Curiouser and curiouser!" cried Alice (she was so much surprised, that for the 
moment she quite forgot how to speak good English). 

- Alice’s Adventures in Wonderland 
Lewis Carroll 
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APPLICATION TO AXIS YM METRIC SUPERSONIC FLOW 


As another illustration of this method, we consider steady, inviscid, supersonic 
flows past axisymmetric bodies. For this purpose consider a family of profiles 
depending on two parameters, thickness and taper. As in Chapter 1, we shall 
summarize the methods used in solving such flows [1], [2]. The equations are written 
in characteristic form as follows: 


s 3 = 0 




(2.1) 

(0 + P(f z) ) a = 

sin2a g i 

(a) - - 

tan 0 tan/z r oc 

(2.2) 

27 

tan 0 + tan/z r 

(8 - m ) e = 

(1 - tan8 

tan/z) 

*-&e , r e 

Y a + tan/z — 
A a r 

(2.3) 


Here the coordinates (a, 3) correspond to the streamlines, a = constant, and the C* 

characteristics, 3 = constant (Figure 1). 8 is the flow deflection angle, /z is the Mach 

angle and s is the entropy. P(/z) is the Prandtl function given by 

P(g)= X^ tan' 1 (X^tan/z) - n , X = (7 + 1)/ (7-1) . (2.4) 

As in the two-dimensional case, these equations are exact and are valid in the 

hypersonic flow regime so long as such real gas effects as disassociation and ionization 
can be ignored. The physical coordinates x,r satisfy the relations [1], [2] 

r <x = x a ten ( 0+ M), r g = Xg tan 0 (2.5) 

As with Chapter 1, the transformation to (a,3) coordinates leaves open two 
arbitrary functions; these are fixed so that the shock is along a = 3 and the body is 
positioned along the line a = 0 (Figure 2). Hence, we have a body fit, shock fit 

coordinate system. The boundary conditions at the body are 

x(0,3) = 3, r (0,3) = f (3,1) , 0 (0,3) = tan 1 (f 0 (3), e) . (2.6) 
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The Rankine-Hugoniot conditions govern the jumps in 6, n and s at the shock. 
Written in terms of the shock angle n, they are given by 


1 (M 2 - 1) t an 2 n - 1 

tan8 = 

tann ( 1+ 2±1 m 2 ) + (1+2^ M 2 ) tan 2 n 


(2.7) 


0.2(1 +-w )(1 + - w) 

9 6 6 
sin *\l = 

(l+w) (1+0.2M ) - (1+ -w) (1+ iw) 

6 6 

s = 2.5 In (1+ -w) + 3.5 In (1+^w) - 3.5 In (l+w) 
6 6 

where 

w = M 2 sin 2 h - 1 

The shock angle h is related to the coordinates as follows 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


r a + r B 

tanh = dr/dx I = I (2.11) 

shock x a + X B shock 

We have assumed a perfect gas with constant specific heats and hence that 

p = p? exp [(7- l)s] (2.12) 

It should be noted that this formulation eliminates the difficulty mentioned in the 
Introduction. For by using the (a,/3) - coordinate system, a quantity such as 

9p 

1 ) 

signifies the variation with ± at fixed a and 13. In particular it gives the variation of 
pressure say fixed at the body or at the shock. This makes the integration of the 
differential equations significantly simpler. 
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VARTATTONAL EQUATIONS 


As in Chapter 1, we differentiate the governing equations with respect to the 
parameter of interest, keeping the coordinates a and 3 held fixed. The differentiation 
although straightforward is tedious. If we write 


30 3s 3x 3r 3g 

e -aT' S -9T' x '8T ,R -17" |, -ar 


(2.13) 


and parametrically differentiate (2.1), (2.2), (2.3) and (2.5) we then obtain, 
S3 = 0 

* nvoo ‘jf^L .^1 Oil* s a 


r a se c 2 Stan 2 g 


7 r(tan9 + tan g) 2 

R = 0 


e - 


r a tan 0tang 


r(tan0 + tang) 2 r < tan6 + tan *0 

r x 0 

3 3 + — ^-“sec 2 0 tang 0- (1- tan 0 tang) — e a = 

L x „ j x„ 


P g(M) % + [P MM (M)g e 


XR0„ , r 3 T 

- — — sec z g tan© + sec 2 g ] 7 


a 


3 


+ (1- tan0 tang) — (Xg X a ) + 


tang 


Re 


rgtang 


R 


(2.14) 


(2.15) 


(2.16) 



= X a tan( 0+ g) + X a sec 2 (0 + g)(0 + 7) 

(2.17) 


£ 

II 

1 

CD 

+ 

$ 

8 

CD 

(2.18) 

At the 

shock the parametrically differentiated equations are 


db 

cosh 

R «( x a + x 3> + R 3( x a + x 3> * X a (r a + rg) 


d T 

( x a + x 3> 2 



- Xe(r a+ r 0 ) 


(2.19) 
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t 


s = 


17.5 3.5 

3.5 

dw 

6 +7w 6+w 

1+w 

dc 

0.2A f- + 2 wl 

- 0.2 fl+ -w 

l 3 18 J 


^ 6 • 


( 2 . 20 ) 


Y = 


A 2 sin2g 


3 18 J dw 

dT 


dw d i) •> 7 i 

where JT = M 2 sin 2n ^ ’ A = W>.2M 2 (l + w) - (1+^w) (1+Jw) 


d9 _ dn 
d€ d€ 


COS 2 0f- 

L sin 


(M 2 -l)tan 2 n -1 


(1+ M 2 )+ (l+^M 2 )1an 2 n 


( 2 . 21 ) 


( 2 . 22 ) 


+ [(1 ,2tl M 2 ) * (l^ M 2 )taA] 2 [(M 2 * 1) sec <7(1+ 2±1 M 2) 

+( 1 + M 2 )tan 2I2 l ‘ [(M 2 -l)tan 2 >7-l](l+ M 2 )sec 2 *l] | 


In the actual integration (2.18-2.22) are applied at the shock 

cc = 0 (2.23) 

At the body a = 0 the appropriate equations are 

m af 0( 0 > € ) 7 

X_= o, R = f — , 9 = — cos 2 e (2.24) 

In writing (2.24) we revert to the general case in which many parameters are being 
considered. Now we simultaneously numerically integrate the non-linear system and the 
variational equations. The calculation of the base flow is second order accurate [2]. The 
calculation of the new flow is first order in space, second order in the parameters of 
interest. The calculation of the two flows is interleaved in that after the flow along 
3 = 0 constant is computed by the base code, the parametric code then calculates the 
exact derivatives in order to obtain the variational flow. 
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RESULTS 


In the numerical calculations discussed, we have taken for f a family of shapes 
given by 

r = 2€x(1-x) + 10x€ (x- 1) (x - — ) c (2.25) 

2 

Here we have taken € to be the thickness ratio based on chord, and c as a scaling 
factor for the taper function. 

Figure 3 shows the effect of changing just the thickness (c = 0). Here we have 
plotted the pressure distribution on the upper surface and the body which, for 
aesthetics, has the lower surface plotted as a reflection of the upper surface. Note that 
the method gives good agreement with the exact solution even when the new thickness 
is 50% more than the base thickness. 


Figure 4 shows the effect of changing a combination of thickness, and taper. Here 
we see that, although the body configurations are markedly different, there is very good 
agreement between the parametrically generated pressure distribution and the exact 
pressure distribution for the new body. 

INVERSE CASE 

The method which has been presented also works quite well in the inverse or 
design problem where the pressure on the body is known, but the shape of the body 
shape is to be determined. 

Using the Bernoulli equation and the perfect gas law one may show [1] 


IL = sin 




y-i 

2 


7- 1 

exp(-y (s + lnrp)) 


7-1 2 7-1 

1+ M - exp(— — (s + ln7p)) 

L / 


■fi 


(2.26) 


This when differentiated, yields 
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d (In p) 
d€ 


(7-1) 2 

2y 


( * + ~*~2 M )exp(w) 


[l-t^M 


exp(w 




s = 


sin 2 g 


= ( 7- l) z ( l+^V 2 ) exp(w) 

[l+^— - exp(w)] 
2 


27 


sin z fi 


Where 


(2.27) 


7-1 

w = -y (s+In p) 


At the body, equation (2.6) is no longer valid since we are attempting to 
determine the shape of the body. Instead we must use (2.5) and (2.26). Therefore, 

the parametrically differentiated equations (2.24) must be replaced by (2.17), (2.18) and 
(2.27). The integration may now proceed as in the direct case [2]. The results of the 
variational calculations are presented in Figures 5 and 6. Notice that even for a 10% 
change in the logarithm of the pressure (corresponding to a 20% increase in thickness), 
tihe difference between the exact body shape and the computed shape is less than 1%. 
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FIGURE CAPTIONS 


Fig. 1 

Fig. 2. 
Fig. 3. 

Fig. 4. 

Fig. 5. 

Fig. 6. 


Body, C+ characteristics, streamlines and C characteristics (dashed) in 
physical (x,y) plane, from [1]. 

Body, C+ characteristics and streamlines in (a, 3) plane. 

Pressure distribution on 25% and 30% thick bodies at M = 6 and 
the respective bodies. 

Pressure distribution on untapered, 25% thick body and 0.10 taper, 
30% thick bodies at M = 6 and the respective bodies. 

Inverse case: Pressure distribution on 25% and 30% thick bodies, 

M = 4 along with generated bodies. Dashed body is computed shape. 

Inverse case: Pressure distribution on 25% and 30% thick bodies, 

M = 6 along with generated bodies. Dashed body is computed shape. 


- 28 - 


PHYSICAL PLANE 



X 

Figure 1 


COMPUTATIONAL PLANE 



- 29 - 


























REFERENCES 


[1] Lewis, T.S. and Sirovich, L., "Approximate and Exact Numerical Computation of 

Supersonic Flow Over an Airfoil", J. Fluid 
Mechanics Vol 112, 1981, pp. 265-282. 


[2] Fong, J. and Sirovich, L. 


'Direct and Inverse Problem in Supersonic 
Axisymmetric Flow", AIAA Journal, 4, 5 May 1986. 


Chapter III 

The Jacobi Matrix Method for General Flows 


Here one of the guinea-pigs cheered, 
officers of the court. 


and was immediately suppressed 


by the 


- Alice’s Adventures in Wonderland 
Lewis Carroii 
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DIRECT DIFFERENCING 


The procedure outlined in Chapters 1 and 2 holds in much greater generality than 
we have considered. The Jacobi matrix technique could also be applied to unsteady 
flows and to viscous flows in three dimensions. However, the method as presented so 
far, has one possible drawback which was alluded to earlier - to obtain the Jacobi 
matrix we must analytically differentiate the relevant equations and boundary conditions. 
In this chapter we propose a procedure which will allow for the calculation of the 
Jacobi matrix by the use of differential approximations. The goal is to obtain the 
Jacobi matrix, and hence be able to calculate a range of solutions in parameter space, 
using the results obtained from solving the nonlinear system (LI) at only two distinct 
values of €.. This differential approach will be applied to the case of two dimensional 
supersonic flow considered in Chapter 1 and to two dimensional subsonic potential flow. 

In the Introduction we said that if u° = u(x ; €q) represented a known solution 
of the base flow then any neighboring flow at some fixed point x is approximated by 


u(x;i) * u° + 



G* ' *o> 


(3.1) 


The obvious first order approximation for 



is 



x, fixed 


u(x;e_) - u° 


*0 


x, fixed 


(3.2) 


9u° 

What this says is that to compute we can take the value of u at the location x 

— k 

in the base flow and subtract it from the value of u found from the perturbed flow 
(i. = + A^) at the same location. In practice, this may require interpolation on one 

computational grid. 

This approach requires special attention at a boundary. In our approach both 
material boundaries and possible shocks are taken to be boundaries and both give rise 
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to locations which change with e.. This would certainly be the case if we chose to 
vary the parameters of a body. 

To be more specific, we would like to be able to use the calculation of pressure 
in the base flow in order to compute the pressure on the new body. Thus the 
formulas (1.25), (1.26) are no longer applicable since they apply at a fixed field point. 
Therefore, to correct (1.26) we must include changes in location of the body due to 
changes in e.. In the interests of simplicity we specify a three dimensional body by 

y = f(x,z;e) (3.3) 


A typical quantity, say pressure, at the new body, which we will specify by Xq, is 
related to the old body Xq in the following way 

9p(Xo;io) 9p(Xo;io) 0f 

P(Xo;±) « Pfeo-ifl) + ' A lk + A€ k (3- 4 ) 

where A £ = i - and 

h = (x 0 ,f(Xo,z 0 ;i 0 ), z o) (3.5) 

and 

Xq = (Xo.fyXo.ZQ-.i),^) (3.6) 


Compare (3.4) with equation (1.26). 

Note that we have related Xq to Xq by placing Xq directly above Xq in the x-z 
plane. Other choices are possible and may be more appropriate in certain cases. 

Equation (3.4) in fact gives us the ability to compute the pressure at the new 

dp 


body, but requires knowledge of the differential coefficients 
obtained from 


de. 


They can be 


3p(Xq ’in ) 


dpte^ df 

p(2^;D - pG^) - — ^ 


9yo 


-k 


9 ^k 


^k 


(3.7) 
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to first order, or 


3_€fl 

3p(Xo;io) pQSo^Cl + € 0 )/2) " P(Xq» — 

^ 

a€ k A ^k 

(3.8) 

It should be noted that the differential determination of the differential coefficient 

9p 9p 

requires not only calculation of the different flow fields, but also of — — . 

9e_ 3y 0 

9p 

Therefore, in the numerical calculation it is necessary to compute — — at the body. 

syo 

This we do by interpolation. 

To illustrate these remarks we return to the case treated in Chapter 1. We 

consider two dimensional supersonic flow at thicknesses of 10% and 10.1% to calculate 
dx dy d9 d/r ds 

, , , , . Using equation (3.8) the resulting derivatives were 

de de de d« de 

then used to compute the pressue distribution on a 15% thick airfoil (Figure 1). Note 
that this pressure distributuion compares favorably with that computed from using the 
Jacobi matrix generated by solving the differential equations (Chapter 1, Figure 3). The 
error between the two computations is less then 1%. 

2D SUBSONIC FLOW 

As a second illustration we apply the Jacobi matrix technique to the potential 
equation for two dimensional compressible flow. The potential equation is derived by 

assuming inviscid, irrotational flow and is valid for subsonic flows and for low transonic 
flows when boundary layer effects can be neglected. 

Since we have implemented the Jacobi technique by modifying Jameson’s 
computer code FL036 we will summarize the derivation of the relevant equations and 
their solution [1],[2]. 

Under the assumption of irrotational flow we may introduce a velocity potential 
<p such that 


£ 

T 


9p(x 0 ; 


3 €_q 


2 ) 9f 


9 y 0 


A ^k 
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(3.9) 


u = <*> x v = 4> y 

The potential satisfies the quasilinear equation 


(a 2 - u 2 ) ^ - 2uv <p xy + (a 2 - v 2 ) <f> yy = 0 (3.10) 


where a is the local speed of sound. When given the ratio of specific heats y, the 
stagnation speed of sound ag and the local speed q = v u “l + y 2 the speed of sound is 
determined by 



(3-11) 


We will consider (3.10) for subsonic flows. (But see Figure 9 for a transonic 
case). 

At the body the flow must satisfy the tangentcy conditionn 



(3.12) 


where n is the normal derivative and the Kutta condition - that the tangential velocity 
is bounded at the trailing edge. In the far field the potential approaches the potential 
of a vortex in compressible flow and a uniform stream. The density and pressure are 
determined by relations 

p7~i = Mia 2 (3.13) 

and 

p7 

p = -1 (3.14) 

7M 2 


The coordinate system used for computation is generated by conformally mapping 
the exterior of the airfoil to the interior of the unit circle. The airfoil itself becomes 
the coordinate line r = 1 (Figure 2). 

Since the far field boundary condition must now be applied at r = 0, where the 
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potential becomes infinite, a reduced potential which removes this singularity is 
introduced by 

cos (0 + a) 

G = <t> - + E(0 + a) (3.15) 

Here a is the angle of attack and 27ZE is the circulation. 

If the modulus of the transformation from the physical plane to the circle plane 
is denoted by H then (3.10) becomes 

0 

(a 2 - u 2 )G ee - 2uvrG r0 + (a 2 - v 2 )r — (rG) 

- 2uv(G e - E) + (u 2 - v 2 )rG r + (u 2 + v 2 ) ( j H 0 + vH f ) = 0 (3.16) 

The u and v are the velocity components in the 9 and r directions, respectively and 
are given by 

r(G 0 - E) - sin(0 + a) r 2 G f - cos (9 + oc) 

u = , v = (3.17) 

H H 

The Neumann boundary condition (3.12) becomes 

G = cos(0 + a) at r = 1 (3.18) 

while the far field condition is 

G = E{0 + a. - tan' 1 [v{ . tan(0 + a) ] } at r = 0 (3.19) 

The circulation is determined by the Kutta condition which requires that the velocity be 
finite at the trailing edge of the airfoil. Here we have H = 0 and <f> 0 = 0 so (3.15) 
reduces to 

E = G 0 - sina at r = 1, 0 = 0 (3.20) 

The details of the calculation of H and of the multigrid solution of (3.16)-3.20) 


- 40 - 


are not essential for our purposes and are discussed in references [1], [3], [4], [5]. The 
important point is that the transformation to the circle plane is conformal so that every 
airfoil in the physical plane is mapped to a circle and every physical flow is mapped to 
the interior of the circle. 

CALCULATION OF THE JACOBI MATRIX 

The equations for computing the Jacobi matrix by finite differences were given by 

equations (3.4)-(3.8). In the subsonic case the only parameter changed was the 

airfoil thickness based on chord. Due to the construction of equation (3.16) the 

quantities which are of interest are the reduced potential G, and the metric H. To 

define the locations x in equation (3.2) we note that in the circle plane the points are 

spaced angularly (0) as 2n/(the number of grid points about airfoil) and radially (r) as 

l/(the number of grid points from airfoil to far field). Therefore, it is natural to 

define the location x by the intersection of these lines. 

The variational flow was computed using essentially the same procedure which was 

used to calculate the pressure in the two dimensional supersonic flow case. In the 

transformed plane we first compute the flow about an airfoil of thickness e 1 and save 

the converged values of G and H. Next we compute the flow about an airfoil of 

thickness € Q . We use these computed values of G and H along with those from the 

dG dH 

run at tnickness €, to compute and using (3.7). G and H for the 

1 de de 

variational flow at thickness e is computed using equation (3.4). 

RESULTS 

Figure 3 shows the results of a parametric calculation using a base airfoil of 
10% thickness based on chord and a second airfoil of 10.1% thickness to predict the 
pressure distribution on a 14% thick airfoil. It should be noted that there is very 
close agreement between the parametric calculation and the solution given by FL036. 

Figure 4 uses a 10% thick airfoil and 10.1% thick airfoil to calculate the flow 
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over 15% thick profile - that is, a profile which is 50% more than the base airfoil. 
Again the agreement is quite good. Figures 5 through 8 show the same calculations 
for flows at different Mach numbers. All show close agreement between the 
parametrically generated solutions and those given by FL036. 

The method breaks down when there is a drastic change in the behavior of the 
solution in the parameter space. This is illustrated in Figure 9. Here the flows about 
the 10% and 10.1% thick airfoils are subsonic but the flow about the 15% thick airfoil 
is supercritical. The method is unable to account for the shock. 
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FIGURE CAPTIONS 


Fig. 1. 

Fig. 2. 
Fig. 3. 

Fig. 4. 

Fig. 5. 

Fig. 6. 

Fig. 7. 

Fig. 8. 

Fig. 9. 


Pressure distribution on 10% and 15% thick airfoils at M = 4 
calculated by direct differencing along with respective bodies. 

Computational plane, from [2]. 

Pressure distribution on 14% thick airfoil, M = 0.75. Solid curve is 
FL036 result, dashed is parametric. Base and new airfoils are also 
shown. 

Pressure distribution on 15% thick airfoil, M = 0.75. Solid curve is 
FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 14% thick airfoil, M = 0.60. Solid curve is 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 15% thich airfoil, M = 0.60. Solid curve if 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 14% thick airfoil, M = 0.45. Solid curve is 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 15% thick airfoil, M = 0.45. Solic curve is 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 15% thick airfoil, M = 0.80. Solid curve is 

FL036 results, dashed is parametric. Base and new airfoils are also 
shown. 
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